"Kohn-Shamification" of the classical density-functional theory of inhomogeneous 
polar molecular liquids with application to liquid hydrogen chloride 
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The Gordian knot of density-functional theories for classical molecular liquids remains finding an 
accurate free-energy functional in terms of the densities of the atomic sites of the molecules. Fol- 
lowing Kohn and Sham, we show how to solve this problem by considering noninteracting molecules 
in a set of effective potentials. This shift in perspective leads to an accurate and computationally 
tractable description in terms of simple three-dimensional functions. We also treat both the linear- 
and saturation- dielectric responses of polar systems, presenting liquid hydrogen chloride as a case 
study. 

PACS numbers: 71.15.Mb, 05.20.Jj, 77.84.Nh 



Introduction — The importance of inhomogeneous po- 
lar molecular liquids, above all water, in the physical 
sciences can hardly be overstated: as solvents they are 
ubiquitous in soft condensed matter physics, biophysics, 
nanophysics, and chemical physics. The associated phe- 
nomena include hydrophobic interactions P, 0] , protein 
folding 0, 01 , the behavior of colloidal suspensions Q, 
and phase transitions of confined liquids [6|, |7| . Due to 
the complex interplay of hydrogen bonding, long-range 
polar interactions, and short-range excluded volume ef- 
fects, developing a tractable physical theory to describe 
the solvent in these systems remains a challenge Q. 

Despite the importance, inherent interest, and exten- 
sive experimental study of polar liquids, most existing 
theories of the inhomogeneous molecular liquid either 
lack the accuracy to describe the above phenomena in 
a quantitatively satisfying way or become computation- 
ally prohibitive when applied to such complex problems. 
The description of liquids from first principles via ab ini- 
tio methods may be accurate, but can only be applied 
to relatively small systems jsij. Even classical molecular 
dynamics requires rather long simulation times to sample 
phase-space sufficiently to extract meaningful thermody- 
namic averages [l^. This latter approach also suffers 
from the well-recognized difficulty of designing potentials 
to describe hydrogen-bonded liquids accurately [U | . 

An alternative approach is to start with the quantum 
mechanical free-energy functional for both electrons and 
nuclei and, by "integrating out" the electrons, to con- 
struct a density functional in terms of atomic site densi- 
ties alone. These "classical" density-functional theories, 
which have been successfully applied to the study of sim- 
ple liquids [13, provide a description of inhomoge- 
neous physical systems that is founded on a number of 
exact theorems [M, 15|- To apply this approach to the 
study of simple liquids, a hard sphere reference system 
is usually augmented by terms that capture weak long- 
range attractive forces [16|. Unfortunately, for most liq- 
uids of interest, a hard sphere reference system is a poor 
starting point because of the strong anisotropic short- 



range interactions arising from the molecular structure 
and effects such as hydrogen bonding. 

To remedy this. Chandler and coworkers 17, [3, [iSl 
introduced a density-functional theory for molecular liq- 
uids in terms of a set of densities, one for each "inter- 
action site" on the molecule (typically atomic centers). 
However, the construction of accurate free-energy func- 
tional in such theories is challenging due to the "inver- 
sion problem" , the difficulty of using only atomic site den- 
sities to express the entropy associated with the geomet- 
ric structure of the molecules. Below, we show how this 
inversion problem can be overcome by a Kohn-Sham-like 
change of variables from site densities to effective poten- 
tials and how the resulting functionals are both compu- 
tationally tractable and can capture the basic underlying 
physics of molecular liquids, including dielectric screen- 
ing effects. 

Kohn-Shamification — The grand free energy f2("*) of 
a noninteracting gas of molecules is well-known as a func- 
tional of the relative potentials ■0Q(r) = (pai^) — t^cn where 
0Q(r) is the site-dependent external potential and /i^ is 
a site-specific chemical potential. 



(1) 



Here, is the reference density at vanishing chemical 
potentials, M is the number of interaction sites on the 
molecule, and s{{ra}), which describes the geometry of 
the molecule, is the intra-molecular distribution function. 
For rigid molecules (zero internal energy) , the entropy of 
the atomic site densities is easily extracted from fi^"*'. 
Thus, to construct an exact density-functional for the 
entropy, even for the interacting system, one only need 
express Jl^"*) as a functional of the site densities. More 
generally, fi^"*) is a key part of the interacting functional 
directly analogous to the noninteracting kinetic energy 
functional Tg [n] of Kohn and Sham . Unfortunately, 
(H]) only gives 0^"'^ as a functional of the t/^a- To ex- 
press ri^'") in terms of site densities requires solution of 
(r) = (5Sl("*^/(5V'Q(r) as a set of M coupled integral 
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equations for the V'a(r) in terms of the n'^^\r). This 
constitutes the "inversion problem" described above. 

Chandler and coworkers have solved this inversion 
problem analytically in terms of an infinite continued 
fraction of convolutions for the rfzatomic case only [3]. 
The unwieldiness of this formal solution and its limitation 
to diatomics, however, has led actual calculations to be 
performed in the united-atom approximation, where, for 
the noninteracting part of the functional, all of the sites 
of a molecule are assumed to coincide at a single point 
figt - Despite its crudeness, the resulting theory predicts 
quite well the correct stable ice structure at nearly the 
correct density [l^. However, by uniting the sites, this 
approximation cannot capture the dielectric response of 
an ideal gas of polar molecules and so provides a poor 
starting point for the study of dielectric effects. 

The key observation which allows us to overcome the 
inversion problem is that the free energy of a noninter- 
acting molecular system, being a very complicated func- 
tional of the site densities, is a very simple functional 
of the relative potentials. This parallels the situation 
in electronic density-functional theory, where there is no 
known accurate functional for Ts[n], the noninteracting 
(kinetic) energy, in terms of the density, but where this 
energy is easily written exactly in terms of single-particle 
orbitals. It was the change in variables from the density 
to the single-particle orbitals of fictitious noninteracting 
particles, as proposed by Kohn and Sham ^d\, which al- 
lowed for the construction of accurate density functionals 
for the interacting electron gas. In the present case, the 
simplicity of evaluation of the grand free energy of nonin- 
teracting molecules as a functional of the relative poten- 
tials suggests an analogous change of variables, now from 
the site densities to a set of effective relative potentials 
in which fictitious noninteracting molecules move. 

Mathematically, a pair of Legendre transformations 
achieves this change of variables. Thermodynamically, 
the Legendre transform fi^"*) ~ / dripaUa equals the 
intrinsic Helmholtz free energy of the noninteracting sys- 
tem. Adding the internal energy U due to inter- molecular 
interactions and the contribution of the physical exter- 
nal potentials (j)a yields the full Helmholtz free energy 
of the interacting system. Finally, the second Legendre 
transformation subtracts J dr nafia to form the full 
interacting grand free energy as a functional of the ipa , 

M „ 

r! = l](™)-^ / d3r(V„(r)-0„(r)-}-^„)rv(r) + C/[n], 

Q = l 

(2) 

where n — {ni(r), riAf (r)}, the set of densities, is 
explicitly a functional of the effective relative poten- 
tials * = {V'i(r-),...,VM(r-)} via n„(r) = ni™'[*](r) = 
jQ(ni) I ^ rpj^g effective relative potentials that min- 
imize then determine the equilibrium site densities 
and allow for the calculation of the various equilibrium 
properties of the liquid. 



Construction of approximate functionals — As usual 
with density-functional theories, the construction of the 
internal energy U is difficult. Again, we follow the lead 
of Kohn and Sham and construct a free-energy func- 
tional that reproduces established results for the homoge- 
neous phase in the limit of vanishing external fields. The 
Ornstein-Zernike (OZ) equation gives information about 
the analytic structure of U, in particular it gives d^U, 
the Hessian of U with respect to the densities, as the 
difference between the inverses of the full and noninter- 



acting correlation function matrices 21|. In the special 
case of the homogeneous phase, translational invariance 
then turns the OZ equation into a simple matrix equation 
in Fourier space, with one component for each field. 

From now on, we will limit the discussion to a spe- 
cial class of liquids, which includes all diatomics and the 
liquid of most interest, water. The liquids in this class 
have the property that the lowest-order term in the long- 
wavelength expansion of the Hessian of U[n] has the form 



Am) 
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where are the partial site charges, e is the macroscopic 
dielectric constant and e^'"' is the dielectric constant of 
a system with intra-molecular correlations only. This re- 
sult may be derived by expanding the interacting and 
noninteracting correlation functions as E + Fk^ + . . . , 
where there are no linear terms in k by rotational in- 
variance, and E and F are matrices of coefficients, with 
one combination of coefficients from F giving the dielec- 
tric constant. In the cases of molecules described either 
by two sites or by three sites with a mirror plane sym- 
metry, one can show that the only combination of the 
coefficients in F which enters the leading-order term in 
the OZ equation corresponds exactly to the bulk dielec- 
tric constant. All liquids composed of such molecules 
thus have the property that the long-wavelength dielec- 
tric response can be built into the free-energy functional 
without any knowledge of the long-wavelength limit of 
the experimental correlation functions (other than e). 

In general, we may view the internal energy U as 
expanded in a power series about the uniform liquid, 
with all higher-order terms gathered into an excess part, 
F*^^, which plays exactly the same role as the exchange- 
correlation functional of Kohn-Sham theory: it ensures 
proper bulk thermodynamic behavior and will ultimately 
be treated in some approximate way. Below, the constant 
term in the power series for U also will be handled as 
part of F'^^ . Finally, the linear term vanishes in the uni- 
form equilibrium state of the system, and the quadratic 
term is given by the Hessian, which ([3]) gives in the long- 



3 



wavelength limit. Putting this all together gives 




C^^{v,v')}n^{r') + F^^[n], (4) 

where Ca-y describes those parts of the Hessian d^U 
which K fails to capture. 

To approximate F'^^ , we first note that, in the case 
of zero external fields, all densities are equal and the 
first quadratic term (K) in ^ vanishes because of charge 
neutrality. Anticipating that the matrix function C will 
be constructed to vanish in the long- wavelength limit, 
F"^^ then captures all the internal energy of the uniform 
phase and can be expressed as = l//^^(n), with V 
being the volume and n the average molecular density. 
Due to the presence of multiple density fields, general- 
izing this expression to the inhomogeneous case is more 
difficult than for the analogous exchange-correlation en- 
ergy in electronic structure theory. Also, because of the 
strong correlations induced by excluded volume effects, 
purely local excess functional fail to describe the liquid 
state [i^]. We therefore approximate F*^^ with a sim- 
plified ansatz in the spirit of weighted density-functional 
theory p^ . but generalized to multiple species, 

F-[n] ^ J d^Y^p^r (E^;%W^ , (5) 

where we introduce the weighted densities n^i^^) = 
/ d^r'(7rro)~^/^ exp(-|r - r'|/rQ)n^(r'), with being a 
parameter ultimately fit to the experimental surface ten- 
sion. To reduce to the correct form in the uniform case 
Pi and h\ must fulfill '}2^P^ ^ ^ a-^id Y,^'=i b\ ^l- 

To capture the behavior of the scalar function f^'^in), 
we use a polynomial fit to various bulk thermodynamic 
conditions. The condition that C vanishes in the long- 
wavelength limit subsequently fixes pi and bl^ and ensures 
that C will be bandwidth limited and thus amenable to 
numerical approximation. For a given rg, this then com- 
pletely specifies our approximation to F'^^. Next, relating 
K+C+d^F'^^ to the density-density correlation functions 
through the OZ relation then gives the matrix function 
C for a given rg. Finally, we can determine tq by adjust- 
ment until calculations of the liquid- vapor interface give 
the correct surface tension. 

Hydrogen chloride — We choose liquid hydrogen chlo- 
ride as a model physical system exhibiting hydrogen 
bonding and for which detailed experimental data are 
available, including site-site correlation functions psi. [23| . 

Assuming rigid intra-molecular bonds, the distribution 
function becomes s(ri,r2) = <5(|ri — — _B)/47ri?^, 
with the gas phase bond length B — 1.275 A taken 
from experiment (23 | and 1 and 2 referring to hydro- 
gen and chlorine, respectively. We then select the par- 
tial charges to yield the experimental gas phase dipole 




distance from capacitor plate [nm] 

FIG. 1: Hydrogen- (thick dotted curves) and chlorine- (thick 
solid curves) site density versus distance from hard wall in 
(a) moderate and (b) high apphed field. Reference results for 
zero field given in both panels (thin solid curves). 

moment [1^, resulting in qi = —q2 — 0.171 e with 
e being the fundamental charge. Then, we approxi- 
mate the simple bulk function f'^^{n) by a fourth-order 
polynomial with coefhcicnts adjusted to (a) reproduce 
the bulk modulus of the liquid phase and (b) to al- 
low for coexistence of stable liquid and vapor phases 
at the appropriate densities, where the latter is treated 
as an ideal gas. Next, for the excess free-energy func- 
tional, the requirement that the matrix C vanishes in 
the long-wavelength limit fixes the integrand in ([5]) to be 
[3r^(ni(r))+4r- ((ni(r) + n2(r))/2)-|-3r-(n2(r))]/10, 
under the assumption that the constant term (in an ex- 
pansion in powers of k^) of the inverse correlation func- 
tion is proportional to that of the noninteracting case. 
Then, we determine C as described above using the par- 
tial structure factor data for the uniform liquid (23| and 
the macroscopic dielectric constant [2^ as experimental 
inputs. Adjusting the smoothing parameter to give the 
experimental surface tension |23] yields rg = 4.55 A. 

Results — As a test case, we study the behavior of 
liquid hydrogen chloride at T=194 K, P=4.5 bar (chosen 
due to availability of experimental correlation functions 
[2^), which is in the liquid part of the phase diagram 
somewhere near the triple point. Figures[2a,b) show the 
equilibrium density profiles which our theory predicts in a 
parallel plate capacitor (described as an infinite square- 
well potential for both species) in both moderate and 
strong applied fields, respectively. For comparison, the 
figure also shows the zero- field density profiles. 

The zero-field profiles exhibit an extended gas phase 
region up to a distance of 1 nm from the plates. This 
"lingering" gas phase region exists because molecules can 
minimize the influence of the repulsive walls by leaving 
the system and entering the reservoir at very little free- 
energy cost. As soon as a relatively weak external elec- 
tric displacement D is applied, however, it becomes fa- 
vorable for the dipolar molecules to enter the capacitor. 
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FIG. 2: Electric saturation fraction (P/Psat) versus ap- 
plied electric displacement (D): density-functional results 
(crosses), non-linear electrostatics (solid line), linear response 
(dashed line). 



change of variables yields a numerically efficient and ac- 
curate density- functional theory for molecular liquids. 
The resulting theory has the computational cost of the 
problem of a noninteracting gas of molecules in a self- 
consistent external potential. We have shown how to con- 
struct a functional in this approach which captures the 
coexistence of liquid and vapor phases, the surface ten- 
sion between these phases, and, for the liquid, the bulk 
mechanical properties, site-site correlation functions, lin- 
ear dielectric response, and self-consistent dielectric sat- 
uration effects. This approach can now be applied to 
liquids in complex environments which can be described 
as external potentials, including nanopores, hydropho- 
bic/hydrophilic molecules, and large biosolutes. Without 
modification, the current approach is ready to be devel- 
oped for the liquid of greatest scientific interest: water. 



largely destroying the gas phase region and resulting in 
an almost rigid shift of the density profiles towards the 
capacitor walls for a large range of fields until the non- 
linear dielectric response regime is reached. Figure [U^a) 
shows an example of this behavior for a relatively large 
field but where the dielectric response is still fairly lin- 
ear (a = /3p-D = 2.68, where p is the molecular dipole 
moment). A dramatic qualitative change occurs in the 
strong field case. Figure mb), computed at a = 7.38, 
shows the typical behavior in the high-field case. At 
large fields, the profile for each species exhibits a sharp 
peak followed by strong oscillations, where the peaks and 
oscillations are separated by one molecular bond length 
(~ 0.12 nm) indicating that they are the result of strong 
orientational ordering of the molecular dipoles. The os- 
cillations of the site densities suggest an induced layering 
of molecules close to the surface, which results from the 
reorganization of the hydrogen-bond network in response 
to the increasing dipolar alignment of the molecules. The 
wavelength of the observed oscillations is about ^ 0.3 nm, 
corresponding to the length of the hydrogen bond . 

The equilibrium site densities also determine the in- 
duced polarization P. Figure [2] compares the polariza- 
tion, computed from our density- functional theory, to the 
result of self-consistently screened nonlinear electrostat- 
ics, computed by solving for the polarization P such that 
P = p(™)(D - aMP), where P(™)(£') is the response 
of a gas of noninteracting dipoles in a local field E and 

= e/(e — 1) — e'^"*-' /(e^^"*) — 1) ensures that the correct 
linear response is recovered when D is weak. Figure [2] 
shows that our density-functional theory not only repro- 
duces the linear regime but also captures self-consistent 
saturation effects. We remark that such results cannot 
be obtained within the united-atom approximation or 
any other approximation that does not take the intra- 
molecular bonding geometry into account. 

In conclusion, we have shown how a Kohn-Sham-like 
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